
% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global Ex Ey Ez Gxy Gxz Gyz vxy vxz vyz
S11=1/Ex;
S22=1/Ey;
S33=1/Ez;
S12=-vxy/Ex;
S13=-vxz/Ex;
S23=-vyz/Ey;
S44=1/Gyz;
S55=1/Gxz;
S66=1/Gxy;
S1=[S11,S12,S13,0,0,0];
S2=[S12,S22,S23,0,0,0];
S3=[S13,S23,S33,0,0,0];
S4=[0,0,0,S44,0,0];
S5=[0,0,0,0,S55,0];
S6=[0,0,0,0,0,S66];
S=[S1;S2;S3;S4;S5;S6];
CC=inv(S);
for ii=0:90
angle=ii;
m=cosd(angle);
n=sind(angle);
T1=[m^2,n^2,0,0,0,2*m*n];
T2=[n^2,m^2,0,0,0,-2*m*n];
T3=[0,0,1,0,0,0];
T4=[0,0,0,m,-n,0];
T5=[0,0,0,n,m,0];
T6=[-m*n,m*n,0,0,0,m^2-n^2];
T=[T1;T2;T3;T4;T5;T6];
TE1=[m^2,n^2,0,0,0,m*n];
TE2=[n^2,m^2,0,0,0,-m*n];
TE3=[0,0,1,0,0,0];
TE4=[0,0,0,m,-n,0];
TE5=[0,0,0,n,m,0];
TE6=[-2*m*n,2*m*n,0,0,0,m^2-n^2];
TE=[TE1;TE2;TE3;TE4;TE5;TE6];
C1=inv(T)*CC*TE;
C=C1;
C=C*(10^9);
pi=4*atan(1); 
rou=1500;
c11=C(1,1);
c12=C(1,2);
c13=C(1,3);
c14=C(1,4);
c15=C(1,5);
c16=C(1,6);
c21=C(2,1);
c22=C(2,2);
c23=C(2,3);
c24=C(2,4);
c25=C(2,5);
c26=C(2,6);
c31=C(3,1);
c32=C(3,2);
c33=C(3,3);
c34=C(3,4);
c35=C(3,5);
c36=C(3,6);
c41=C(4,1);
c42=C(4,2);
c43=C(4,3);
c44=C(4,4);
c45=C(4,5);
c46=C(4,6);
c51=C(5,1);
c52=C(5,2);
c53=C(5,3);
c54=C(5,4);
c55=C(5,5);
c56=C(5,6);
c61=C(6,1);
c62=C(6,2);
c63=C(6,3);
c64=C(6,4);
c65=C(6,5);
c66=C(6,6);
N=360;
pahse_angle=zeros(1,N);
vs1=zeros(1,N);
vs2=zeros(1,N);
vp=zeros(1,N); 
g_vs1=zeros(1,N);
g_vs2=zeros(1,N);
g_vp=zeros(1,N);
g_angle_vs1=zeros(1,N);
g_angle_vs2=zeros(1,N);
g_angle_vp=zeros(1,N);
for i=1:1:N
    sita=pi/180*i;
    pahse_angle(i)=sita;
    n2=cos(sita);
    n3=sin(sita);
    n1=0;
    %Christoffel矩阵
    tao11=c11*(n1^2)+c66*(n2^2)+c55*(n3^2)+2*c56*n2*n3+2*c15*n1*n3+2*c16*n1*n2;
    tao22=c66*(n1^2)+c22*(n2^2)+c44*(n3^2)+2*c24*n2*n3+2*c46*n1*n3+2*c26*n1*n2;
    tao33=c55*(n1^2)+c44*(n2^2)+c33*(n3^2)+2*c34*n2*n3+2*c35*n1*n3+2*c45*n1*n2;
    tao12=c16*(n1^2)+c26*(n2^2)+c45*(n3^2)+(c46+c25)*n2*n3+(c14+c56)*n1*n3+(c12+c66)*n1*n2;
    tao13=c15*(n1^2)+c46*(n2^2)+c35*(n3^2)+(c45+c36)*n2*n3+(c13+c55)*n1*n3+(c14+c56)*n1*n2;
    tao23=c56*(n1^2)+c24*(n2^2)+c34*(n3^2)+(c44+c23)*n2*n3+(c36+c45)*n1*n3+(c46+c25)*n1*n2;
    a=[tao11 tao12 tao13;tao12 tao22 tao23;tao13 tao23 tao33];
%     syms v
%     A=[tao11-rou*v,tao12,tao13;tao12,tao22-rou*v,tao23;tao13,tao23,tao33-rou*v];
%     D=det(A);
%     f=factor(D);
%     V=solve(D);
%     V=sqrt(double(V));
    %求特征值
    %[V,D] = eig(A) 返回特征值的对角矩阵 D 和矩阵 V，其列是对应的右特征向量，使得 A*V = V*D。
    [x,y]=eig(a);
    %特征值
    vs1(1,i)=sqrt(y(1,1)/rou);
    vs2(1,i)=sqrt(y(2,2)/rou);
    vp(1,i)=sqrt(y(3,3)/rou);
    %特征值对应的矩阵
    p_vs1(1,i)=x(1,1);p_vs1(2,i)=x(2,1);p_vs1(3,i)=x(3,1);
    p_vs2(1,i)=x(1,2);p_vs2(2,i)=x(2,2);p_vs2(3,i)=x(3,2);
    p_vp(1,i)=x(1,3);p_vp(2,i)=x(2,3);p_vp(3,i)=x(3,3);

for i=1:1:N-1
    g_vs1(1,i)=sqrt(vs1(1,i)^2+(vs1(1,i+1)-vs1(1,i))^2/(pi/180)^2);
    g_vs2(1,i)=sqrt(vs2(1,i)^2+(vs2(1,i+1)-vs2(1,i))^2/(pi/180)^2);
    g_vp(1,i)=sqrt(vp(1,i)^2+(vp(1,i+1)-vp(1,i))^2/(pi/180)^2);

    g_angle_vs1(1,i)=atan(1/vs1(1,i)*(vs1(1,i+1)-vs1(1,i))/(pi/180))+pi/180*i;
    g_angle_vs2(1,i)=atan(1/vs2(1,i)*(vs2(1,i+1)-vs2(1,i))/(pi/180))+pi/180*i;
    g_angle_vp(1,i)=atan(1/vp(1,i)*(vp(1,i+1)-vp(1,i))/(pi/180))+pi/180*i;
end
g_vs1(1,N)=vs1(1,N);
g_vs2(1,N)=vs2(1,N);
g_vp(1,N)=vp(1,N);
g_angle_vs1(1,N)=pi/180*N;
g_angle_vs2(1,N)=pi/180*N;
g_angle_vp(1,N)=pi/180*N;
end
% figure
% polar(pahse_angle,vp,'-b');
% hold on;
% polar(pahse_angle,vs1,'--r');
% hold on;
% polar(pahse_angle,vs2,'--g');
% hold off
vout(ii+1,:)=vp;
end


end